#include "ParticleSource2D_Load.h"
#include "../Parallel/ParallelMPI2D.h"
#include "../Field/Field2D.h"

#include <cmath>
#include <iomanip>
#include <algorithm>
#include <ostream>
#include <sstream>

using namespace std;


//Constructor
ParticleSource2D_Load::ParticleSource2D_Load(PicParams* params_, ParallelMPI* pmpi_, PhysicalField* fields_, vector<Species*>* vecSpecies_, ParticleSourceParameter ps_param_):
    ParticleSource2D(params_, pmpi_, fields_, vecSpecies_, ps_param_)
{
    dt_ov_dx = params->dt / params->cell_length[0];
    dt = params->dt;

    YZArea = 1.0;

    //define the x-direction range of source region
    //the MPI domain is not in the source region
    if(pmpi->getDomainLocalMax(0) <= ps_param_.load_position_start[0] || pmpi->getDomainLocalMin(0) >= ps_param_.load_position_end[0]) {
        load_cell_start[0] = 0;
        load_cell_end[0] = 0;

        load_cell_start[1] = 0;
        load_cell_end[1] = 0;
    }
    //the left end of source region is in the MPI domain
    else if(pmpi->getDomainLocalMin(0) <= ps_param_.load_position_start[0] && pmpi->getDomainLocalMax(0) > ps_param_.load_position_start[0]
    && pmpi->getDomainLocalMax(0) <= ps_param_.load_position_end[0]) {
        double temp = (ps_param_.load_position_start[0] - pmpi->getDomainLocalMin(0)) / params->cell_length[0];
        load_cell_start[0] = temp + 0.001;
        if(load_cell_start[0] < 0) { load_cell_start[0] = 0; }
        //number_particle_in_each_cell[load_cell_start[0]] /=
        //( params->cell_length[0] / (( pmpi->getDomainLocalMin(0) + (load_cell_start[0]+1)*params->cell_length[0] ) - ps_param_.load_position_start[0]) );
        load_cell_end[0] = params->dims_space[0] - 1;
    }
    //the right end of source region is in the MPI domain
    else if(pmpi->getDomainLocalMin(0) < ps_param_.load_position_end[0] && pmpi->getDomainLocalMax(0) >= ps_param_.load_position_end[0]
    && pmpi->getDomainLocalMin(0) >= ps_param_.load_position_start[0]) {
        load_cell_start[0] = 0;
        load_cell_end[0] = (ps_param_.load_position_end[0] - pmpi->getDomainLocalMin(0)) / params->cell_length[0];
        if(load_cell_end[0] > params->dims_space[0] - 1) { load_cell_end[0] = params->dims_space[0] - 1; }
        //number_particle_in_each_cell[load_cell_end[0]] /=
        //( params->cell_length[0] / (ps_param_.load_position_end[0] - ( pmpi->getDomainLocalMin(0) + load_cell_end[0]*params->cell_length[0] )) );
    }
    //the whole MPI domain is in the source region
    else if(pmpi->getDomainLocalMin(0) >= ps_param_.load_position_start[0] && pmpi->getDomainLocalMax(0) <= ps_param_.load_position_end[0]) {
        load_cell_start[0] = 0;
        load_cell_end[0] = params->dims_space[0] - 1;
    }
    //the whole MPI domain source region is in the bin
    else if(pmpi->getDomainLocalMin(0) < ps_param_.load_position_start[0] && pmpi->getDomainLocalMax(0) > ps_param_.load_position_end[0]) {
        //load_cell_start[0] = (ps_param_.load_position_start[0] - pmpi->getDomainLocalMin(0)) / params->cell_length[0];
        //number_particle_in_each_cell[load_cell_start[0]] /=
        //( params->cell_length[0] / (( pmpi->getDomainLocalMin(0) + (load_cell_start[0]+1)*params->cell_length[0] ) - ps_param_.load_position_start[0]) );
        //load_cell_end[0] = (ps_param_.load_position_end[0] - pmpi->getDomainLocalMin(0)) / params->cell_length[0];
        //number_particle_in_each_cell[load_cell_end[0]] /=
        //( params->cell_length[0] / (ps_param_.load_position_end[0] - ( pmpi->getDomainLocalMin(0) + load_cell_start[0]*params->cell_length[0] )) );
        double temp = (ps_param_.load_position_start[0] - pmpi->getDomainLocalMin(0)) / params->cell_length[0];
        load_cell_start[0] = temp + 0.001;
        if(load_cell_start[0] < 0) { load_cell_start[0] = 0; }
        load_cell_end[0] = (ps_param_.load_position_end[0] - pmpi->getDomainLocalMin(0)) / params->cell_length[0];
        if(load_cell_end[0] > params->dims_space[0] - 1) { load_cell_end[0] = params->dims_space[0] - 1; }

    }

    //define the y-direction range of source region
    //the MPI domain is not in the source region
    if(pmpi->getDomainLocalMax(1) <= ps_param_.load_position_start[1] || pmpi->getDomainLocalMin(1) >= ps_param_.load_position_end[1]) {
        load_cell_start[0] = 0;
        load_cell_end[0] = 0;

        load_cell_start[1] = 0;
        load_cell_end[1] = 0;
    }
    //the down end of source region is in the MPI domain
    else if(pmpi->getDomainLocalMin(1) <= ps_param_.load_position_start[1] && pmpi->getDomainLocalMax(1) > ps_param_.load_position_start[1]
    && pmpi->getDomainLocalMax(1) <= ps_param_.load_position_end[1]) {
        double temp = (ps_param_.load_position_start[1] - pmpi->getDomainLocalMin(1)) / params->cell_length[1];
        load_cell_start[1] = temp + 0.001;
        if(load_cell_start[1] < 0) { load_cell_start[1] = 0; }
        //number_particle_in_each_cell[load_cell_start[0]] /=
        //( params->cell_length[0] / (( pmpi->getDomainLocalMin(0) + (load_cell_start[0]+1)*params->cell_length[0] ) - ps_param_.load_position_start[0]) );
        load_cell_end[1] = params->dims_space[1] - 1;
    }
    //the up end of source region is in the MPI domain
    else if(pmpi->getDomainLocalMin(1) < ps_param_.load_position_end[1] && pmpi->getDomainLocalMax(1) >= ps_param_.load_position_end[1]
    && pmpi->getDomainLocalMin(1) >= ps_param_.load_position_start[1]) {
        load_cell_start[1] = 0;
        load_cell_end[1] = (ps_param_.load_position_end[1] - pmpi->getDomainLocalMin(1)) / params->cell_length[1];
        if(load_cell_end[1] > params->dims_space[1] - 1) { load_cell_end[1] = params->dims_space[1] - 1; }
        //number_particle_in_each_cell[load_cell_end[0]] /=
        //( params->cell_length[0] / (ps_param_.load_position_end[0] - ( pmpi->getDomainLocalMin(0) + load_cell_end[0]*params->cell_length[0] )) );
    }
    //the whole MPI domain is in the source region
    else if(pmpi->getDomainLocalMin(1) >= ps_param_.load_position_start[1] && pmpi->getDomainLocalMax(1) <= ps_param_.load_position_end[1]) {
        load_cell_start[1] = 0;
        load_cell_end[1] = params->dims_space[1] - 1;
    }
    //the whole MPI domain source region is in the bin
    else if(pmpi->getDomainLocalMin(1) < ps_param_.load_position_start[1] && pmpi->getDomainLocalMax(1) > ps_param_.load_position_end[1]) {
        //load_cell_start[0] = (ps_param_.load_position_start[0] - pmpi->getDomainLocalMin(0)) / params->cell_length[0];
        //number_particle_in_each_cell[load_cell_start[0]] /=
        //( params->cell_length[0] / (( pmpi->getDomainLocalMin(0) + (load_cell_start[0]+1)*params->cell_length[0] ) - ps_param_.load_position_start[0]) );
        //load_cell_end[0] = (ps_param_.load_position_end[0] - pmpi->getDomainLocalMin(0)) / params->cell_length[0];
        //number_particle_in_each_cell[load_cell_end[0]] /=
        //( params->cell_length[0] / (ps_param_.load_position_end[0] - ( pmpi->getDomainLocalMin(0) + load_cell_start[0]*params->cell_length[0] )) );
        double temp = (ps_param_.load_position_start[1] - pmpi->getDomainLocalMin(1)) / params->cell_length[1];
        load_cell_start[1] = temp + 0.001;
        if(load_cell_start[1] < 0) { load_cell_start[1] = 0; }
        load_cell_end[1] = (ps_param_.load_position_end[1] - pmpi->getDomainLocalMin(1)) / params->cell_length[1];
        if(load_cell_end[1] > params->dims_space[1] - 1) { load_cell_end[1] = params->dims_space[1] - 1; }

    }

    //load particles by number_density and Temperature
    numPart_in_each_cell = loadDensity / params->species_param[species1].weight;
    if(loadKind == "nT") {
        for(int ibin = 0; ibin < number_particle_in_each_cell.size(); ibin++)
        {
            number_particle_in_each_cell[ibin] = (load_cell_end[1] - load_cell_start[1] + 1) * numPart_in_each_cell;
            //cout<<"number_particle_in_each_cell "<<number_particle_in_each_cell[ibin]<<endl;
        }
    }
    else if(loadKind == "dn")
    {
        load_step = 1.0 + load_number * params->species_param[species1].weight / (load_dn * params->dt);
        load_number_rem = load_dn * load_step * params->dt / params->species_param[species1].weight - load_number;
        load_number_rem_tot = 0.0;
        log_info<<"load_step = "<<load_step;
    }


}

ParticleSource2D_Load::~ParticleSource2D_Load()
{

}


void ParticleSource2D_Load::emitLoad(int itime)
{

}


/* Calculates the ParticleSource for a given ParticleSource object
void ParticleSource2D_Load::emitLoad(PicParams* params_, ParallelMPI* pmpi_, vector<Species*>& vecSpecies, int itime, PhysicalField* fields)
{
    Species   *s1;
    Particles *p1;
    int iPart, nPart;
    vector <double> cell_length;            //cell_length for initialize position, maybe not equal to real cell lenghth
    vector<double> max_jutt_cumul;
    double *indexes=new double[params->n_dim_field];
    double *temp=new double[3];
    double *vel=new double[3];

    if(everyTime == 0 && itime > 1) { return; }
    if(itime % step_update != 0){ return; }
    if(loadKind == "nT" && load_cell_end[0] != load_cell_start[0] && load_cell_end[1] != load_cell_start[1])
    {
        s1 = vecSpecies[species1];
        p1 = &(s1->particles);

        cell_length.resize(params->n_dim_field);
        max_jutt_cumul.resize(0);
        temp[0] = load_temperature;
        temp[1] = load_temperature;
        temp[2] = load_temperature;
        vel[0] = mean_velocity[0];
        vel[1] = mean_velocity[1];
        vel[2] = mean_velocity[2];

        indexes_of_particles_to_erase.clear();
        for(int ibin = 0; ibin < count_of_particles_to_insert.size(); ibin++ )
        {
            count_of_particles_to_insert[ibin] = 0;
        }

        //erase unnecessary paritcles in source region
        for(int ibin=load_cell_start[0]; ibin<=load_cell_end[0]; ibin++)
        {
            for(int iPart = s1->bmin[ibin]; iPart < s1->bmax[ibin]; iPart++)
            {
                if( p1->position(1,iPart) > ps_param_.load_position_start[1] && p1->position(1,iPart) < ps_param_.load_position_end[1] ) {
                    indexes_of_particles_to_erase.push_back(iPart);
                }
            }
            count_of_particles_to_insert[ibin] = number_particle_in_each_cell[ibin];
        }
        s1->erase_particles_from_bins(indexes_of_particles_to_erase);

        //insert lost particles into bins
        new_particles.clear();
        for(int ibin = load_cell_start[0]; ibin <= load_cell_end[0]; ibin++ )
        {
            new_particles.create_particles(count_of_particles_to_insert[ibin]);
        }
        s1->insert_particles_to_bins(new_particles, count_of_particles_to_insert);

        //re-initialize paritcles in source region
        for(int ibin = load_cell_start[0]; ibin <= load_cell_end[0]; ibin++)
        {
            //Because bin is only in x-direction, y-direction has no bin
            //The "Bin" in load_cell_start[1] and load_cell_end[1] is just to be consistent with the x-direction
            for(int j = load_cell_start[1]; j <= load_cell_end[1]; j++)
            {
                iPart = s1->bmax[ibin] - (load_cell_end[1] - j + 1) * numPart_in_each_cell ;
                nPart = numPart_in_each_cell;
                cell_length[0] = params->cell_length[0];
                cell_length[1] = params->cell_length[1];
                indexes[0] = pmpi->getDomainLocalMin(0) + ibin*params->cell_length[0];
                indexes[1] = pmpi->getDomainLocalMin(1) + j   *params->cell_length[1];

                s1->initPosition(nPart, iPart, indexes, params->n_dim_field,
                             cell_length, s1->species_param.init_position_type);

                s1->initMomentum(nPart,iPart, temp, vel,
                             s1->species_param.init_velocity_type, max_jutt_cumul, params);

                s1->initWeight_constant(nPart, species1, iPart, s1->species_param.weight);
                s1->initCharge(nPart, species1, iPart, s1->species_param.charge);
            }
        }
    }
    else if(loadKind == "dn" && itime%load_step == 0 && load_cell_end[0] != load_cell_start[0])
    {
        s1 = vecSpecies[species1];
        p1 = &(s1->particles);

        cell_length.resize(params->n_dim_field);
        max_jutt_cumul.resize(0);
        temp[0] = load_temperature;
        temp[1] = load_temperature;
        temp[2] = load_temperature;
        vel[0] = mean_velocity[0];
        vel[1] = mean_velocity[1];
        vel[2] = mean_velocity[2];

        double load_number_temp;
        load_number_rem_tot += load_number_rem;
        load_number_temp = load_number;
        if(load_number_rem_tot > 1.0)
        {
            load_number_temp = load_number + 1;
            load_number_rem_tot -= 1.0;
        }
        numPart_in_each_cell = load_number_temp;

        for(int ibin = 0; ibin < count_of_particles_to_insert.size(); ibin++ )
        {
            count_of_particles_to_insert[ibin] = 0;
            if(ibin >= load_cell_start[0] && ibin <= load_cell_end[0])
            {
                count_of_particles_to_insert[ibin] = (load_cell_end[1] - load_cell_start[1] + 1) * numPart_in_each_cell;
            }
        }

        new_particles.clear();
        for(int ibin = load_cell_start[0]; ibin <= load_cell_end[0]; ibin++ )
        {
            new_particles.create_particles(count_of_particles_to_insert[ibin]);
        }
        s1->insert_particles_to_bins(new_particles, count_of_particles_to_insert);

        //re-initialize paritcles in source region
        for(int ibin=load_cell_start[0]; ibin<=load_cell_end[0]; ibin++)
        {
            //Because bin is only in x-direction, y-direction has no bin
            //The "Bin" in load_cell_start[1] and load_cell_end[1] is just to be consistent with the x-direction
            for(int j = load_cell_start[1]; j <= load_cell_end[1]; j++)
            {
                iPart = s1->bmax[ibin] - (load_cell_end[1] - j + 1) * numPart_in_each_cell ;
                nPart = numPart_in_each_cell;
                cell_length[0] = params->cell_length[0];
                cell_length[1] = params->cell_length[1];
                indexes[0] = pmpi->getDomainLocalMin(0) + ibin*params->cell_length[0];
                indexes[1] = pmpi->getDomainLocalMin(1) + j   *params->cell_length[1];

                s1->initPosition(nPart, iPart, indexes, params->n_dim_field,
                             cell_length, s1->species_param.init_position_type);

                s1->initMomentum(nPart,iPart, temp, vel,
                             s1->species_param.init_velocity_type, max_jutt_cumul, params);

                s1->initWeight_constant(nPart, species1, iPart, s1->species_param.weight);
                s1->initCharge(nPart, species1, iPart, s1->species_param.charge);
            }
        }
    }

    delete [] indexes;
    delete [] temp;
    delete [] vel;
}
*/
